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Abstract 

In this paper we present an efficient discretization method for the solution of the unsteady incompressible 
Navier-Stokes equations based on a high order (Hybrid) Discontinuous Galerkin formulation. The crucial 
component for the efficiency of the discretization method is the disctinction between stiff linear parts and 
less stiff non-linear parts with respect to their temporal and spatial treatment. 

Exploiting the flexibility of operator-splitting time integration schemes we combine two spatial discretiza¬ 
tions which are tailored for two simpler sub-problems: a corresponding hyperbolic transport problem and 
an unsteady Stokes problem. 

For the hyperbolic transport problem a spatial discretization with an Upwind Discontinuous Galerkin 
method and an explicit treatment in the time integration scheme is rather natural and allows for an efhcient 
implementation. The treatment of the Stokes part involves the solution of linear systems. In this case a dis¬ 
cretization with Hybrid Discontinuous Galerkin methods is better suited. We consider such a discretization 
for the Stokes part with two important features: JJ(div)-conforming finite elements to garantuee exactly 
divergence-free velocity solutions and a projection operator which reduces the number of globally coupled 
unknowns. We present the method, discuss implementational aspects and demonstrate the performance on 
two and three dimensional benchmark problems. 

Keywords: Navier-Stokes equations, Hybrid Discontinuous Galerkin Methods, H(div)-conforming Finite 
Elements, exactly divergence-free, Operator-Splitting, reduced stabilization 


1. Introduction 

1.1. Problem statement 

We consider the numerical solution of the unsteady incompressible Navier Stokes equations in a velocity- 
pressure formulation: 

J A\v{—vVu -\- u®u -\- pi) = / in D , . 

( div u = 0 in D 

with boundary conditions u = ud onT^i C dfl and (izVu — pi) • n = 0 on Tout = 9D \ Td. Here, v = const 
is the kinematic viscosity, u the velocity, p the pressure, and / is an external body force. We consider 
a discretization to Q with high order (Hybrid) Discontinuous Galerkin (DG) finite element methods for 
complex geometries with underlying meshes consisting of possibly curved tetrahedra, hexahedra, prisms and 
pyramids. 


Email addresses: christoph.lehreiifeld@uni-inuenster.de (Christoph Lehrenfeld), joachim.schoeberl@tu-wien.ac.at 
(Joachim Schoberl) 


Preprint submitted to Comp. Meth. Appl. Mech. Eng. 


April 27, 2016 





1.2. Literature 

Since its introduction in the paper of Reed and Hill [1], DG methods have been developed and used 
for hyperbolic problems [a m m El E] and later extended to second-order elliptic (and parabolic) problems 
[71 [HI El uni- DG methods are specifically popular for flow problems mi EH. For incompressible flows 
different finite element methods have been discussed in the literature [niiaiii]- Among these methods 
several DG formulations have been considered, e.g. nnuniiiHiiiiiini- 

DG methods provide a flexibility which can be utilized for different purposes. In our case, the motivation 
to consider this type of discretizations is twofold. First, relevant flow problems that can be modeled with 
the incompressible Navier-Stokes equations are often convection dominated. Using the Upwind mechanism 
DG methods offer a natural way to devise stable discretizations of (dominating) convection. Secondly, 
abandoning -conforming finite element spaces, DG methods allow to consider (only) iJ(div)-conforming 
finite elements (see [HI]) for (Navier-)Stokes problems. This is attractive as it facilitates the design of a 
discretization with exactly solenoidal solution which in turn implies energy-stability of the discretization for 
the Navier-Stokes problem, cf. [T^[2(71 [22] . 

Gompared to Gontinuous Galerkin (GG) methods the number of degrees of freedom of Discontinuous 
Galerkin methods increases significantly. This drawdack is often outwayed by the advantages of the method. 
However, when it comes to solving linear systems Discontinuous Galerkin methods suffer most from drasti¬ 
cally increased globally coupled degrees of freedom. An approach to compensate for this is the concept of 
Hybridization where additional unknowns on element interfaces are introduced. This increases the number 
of degrees of freedom, but introduces two advantages. First, the global couplings are reduced and secondly, 
the structure of the couplings allows to apply static condensation for the element unknowns. The con¬ 
cept has originally been introduced in the context of mixed finite element methods, cf. |21j . In the last 
decade Hybridization has also been applied to Discontinuous Galerkin method for a variety of problems, 
cf. [231 El EH [261127]. We also mention, that recently similar concepts are also known in the literature 
under the name “Hybridized Weak Galerkin” methods [HH] where a slightly different framework is used to 
derive the methods. With the emphasis to treat general polyhedral meshes “Hybrid High-Order” (HHO) 
methods EllSniEI] have recently been introduced where a combination of Hybrid (or hybridized) methods 
and higher order spaces is considered. 

Goncerning HDG discretizations for incompressible flows, we also mention the papers [SHIESIEI]. Fur¬ 
ther, different approaches to implement exact incompressible finite element solutions to the Stokes problem 
using Hybridization have been investigated in [5H [551155] . In m the Stokes-Brinkman problem, the Stokes 
problem with an additional zero order term, has been discretized using a hybridized H(div)-conforming 
finite element formulation. 

An interesting and often praised aspect of Hybrid mixed finite element methods is the fact that a post¬ 
processing step can be used to reconstruct interior unknowns of an increased higher order accuracy for 
elliptic problems m- The same is also possible for HDG methods in mixed formulation which is another 
advantage over conventional DG methods, see e.g. [55]. For elliptic problems an accuracy of order k + 2 
in the volume can be achieved if polynomials of order k are used on the element interfaces. One main 
contribution in this paper is the introduction of a projection operator which achieves the same effect. In 
[22] we already discussed this operator which has recently also been addressed under the name “reduced 
stabilization” in [31 EHlIinillT]. In these papers the construction of the corresponding operator is only 
possible in two dimensions using special integration rules. We explain how the operator can be implemented 
in a more general setting. 

The efficiency and practicability of HDG methods is rarely addressed in the literature. Interesting 
exception are [42| [43] where the computational cost of the method is compared to GG and DG methods. 

The nature of diffusive (viscous) and convective terms appearing in convection-diffusion-type problems 
have a substantially different character. Discretizations of problems involving only one of the two mechanisms 
would typically lead to different spatial and temporal discretizations. It is therefore often desirable to consider 
operator-splitting time discretization schemes which allow for a separate treatment of the different operators 
as in [44l [45] |46l im IlHl [49] . In these papers the spatial discretization for the stiff (diffusion/Stokes) and the 
non-stiff (convection) operators is typically the same. We consider a different treatment of the operators 
with respect to their temporal and spatial discretization. 
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1.3. The concept: Efficiency through operator splitting 

A crucial ingredient in the considered discretization is the fact that the convection and the Stokes 
problem are separated by means of an operator-splitting method. The operator-splitting method is chosen 
such that only operator evaluations of the convection operator are required so that the time integration 
scheme is explicit in terms of the convection operator. This is often affordable as the time step restrictions 
following from the convection operator are the least restrictive ones in the Navier-Stokes problem. Moreover, 
the convection operator is non-linear, s.t. the set up of linear systems of equations and corresponding 
preconditioners or solvers would have to take place every time step which renders implicit approaches for 
the convection very expensive. The Stokes operator is dealt with implicitly, i.e. it appears in linear systems 
in every time step. Due to the differential-algebraic structure of the Navier-Stokes equations this is necessary 
w.r.t. the pressure and the incompressibility constraint. As completely explicit handling of the viscosity 
terms would introduce severe time step restrictions it is also advisable to treat viscous forces implicit. Note 
that the Stokes operator is time-independent such that the setup of linear systems and preconditioners or 
solvers can be done once and re-used in every time step. 

The spatial discretizations for the different operators are designed differently as the treatment of the 
convection term can be optimized for operator evaluations while the discretization of the Stokes operator 
has to provide an efficient handling of implicit solution steps, i.e. the solution of corresponding linear 
systems. This different treatment reflects in the use of two different finite element spaces which are used: 
An i?(div)-conforming Hybrid DG space for the velocity-pressure-pair for solving Stokes problems and a 
DG space for handling convection. 

1 . 4 . Main contributions and structure of this paper 

In this paper we introduce a new discretization method for the incompressible Navier-Stokes equations. 
The discretization is based on a decomposition of the problem into the (unsteady) Stokes problem and 
a hyperbolic transport problem. This decomposition reflects in the use of appropriate operator splitting 
methods in the time discretization and the use of different finite element spaces for the different spatial 
operators. 

While we use a rather standard DG formulation for the spatial discretization of the hyperbolic transport 
problem, we consider a new HDG formulation for the spatial discretization of the Stokes-type problem. We 
use H (div)-conforming functions of element-wise polynomials of degree k for the velocity and discontinuous 
pressure functions of degree k — 1. Due to the introduction of additional unknowns of degree A: — 1 on 
the facets for the approximation of the tangential trace of the velocity static condensation allows to reduce 
the unknowns and their couplings significantly. The resulting globally coupled unknowns correspond to the 
approximation of the normal trace of the velocity on the facets by order k functions, the tangential trace of 
the velocity by order k — 1 functions and the approximation of the pressure field by piecewise constants. 

The main contributions of this paper are: 

1. Introduction of a new iJ(div)-conforming high order accurate HDG method with a projected jumps 
formulation (also known as reduced stabilization or reduced-order HDG) for the solution of Stokes-type 
problems. 

2. Presentation of a combined spatial discretization for the Navier-Stokes equations based on a standard 
Upwind DG formulation for the hyperbolic transport problem with the new HDG method for Stokes- 
Brinkman problems. 

3. Discussion of operator-splitting time integration schemes which restrict solutions of linear systems to 
Stokes-Brinkman problems. 

In section we discuss the discretization of spatial operators. The discussion is divided into two parts, the 
discretization of the Stokes-Brinkman problem, the problem which involves all relevant spatial operators 
except for the convection and the discretization of the convection operator. For the former part we consider 
an HDG discretization, for the latter part a standard DG discretization. In section operator-splitting 
time integration schemes are discussed which are tailored for such a situation. Finally, in section]^ we give 
numerical examples which demonstrate the accuracy and the performance of the method. 
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1.5. Preliminaries 

Before describing the methods, we introduce some basic notation and assumptions: O is an open bounded 
domain in with a Lipschitz boundary F. It is decomposed into a shape regular partition 7^ of consisting 
of elements T which are (curved) simplices, quadrilaterals, hexahedrals, prisms or pyramids. For ease of 
presentation we assume that all elements are not curved and consider only homogeneous Dirichlet boundary 
conditions. The element interfaces and element boundaries coinciding with the domain boundary are called 
facets. The set of those facets F is denoted by and there holds UreTh ~ 

In the sequel we distinguish functions with support only on facets indicated by a subscript F and those 
with support also on the volume elements which is indicated by a subscript T. Compositions of both types 
are used for the HDG discretization of the velocity which is denoted by u = {ut,uf). 

For vector-valued functions the superscripts t denotes the application of the tangential projection: v* = 
V — (v-n) ■ n G The index k which describes the polynomial degree of the finite element approximation 
at many places through out the paper is an arbitrary but fixed positive integer number. 

We identify finite element functions u G Xh with their representation in terms of coefficient vectors 
uGK^^,s.t. u = iLpi for a corresponding basis {(fi} of Xh and Nx = dim(Xft,). A (generic) bilinear 

form Qh ■ Xh x Y/j —)■ M is identified with the matrix G G g 


2. DG/HDG spatial discretization 


In this section we introduce the spatial discretization for the Stokes operator, the convection operator 
and transfer operations between both. First, we introduce the iF(div)-conforming Hybrid DG discretization 
of the Stokes-Brinkman problem, i.e. the stationary Navier-Stokes problem without convection and an 
additional zero order term in section in section |2.1[ This discretization is improved significantly. Further 
on, a modification of the discretization using the idea of projected jumps (also known as reduced stabilization 
or reduced-order HDG) is presented in section 2.2 including an explanation of how the projection operator 


can be realized. For the convection part of the Navier-Stokes problem we consider a DG discretization using 
standard approaches. This is discussed in section 2.3 As the discretization spaces for the Stokes part and 
the convection part are different, we present transfer operations between the spaces in section 2.4 which 
allow us to finally formulate the semi-discrete problem in section |2.5[ 


2.1. F[{div)-conforming HDG formulation of the Stokes-Brinkman problem 

In this part we consider the discretization of the Stokes part of the Navier-Stokes problem. For simplicity 
we restrict the discussion to homogeneous Dirichlet boundary conditions. We present the method and 
elaborate on important properties. For an error analysis of the method we refer to [55]. The reaction 
term corresponding to an inertia term stemming from an implicit time integration scheme is further 
incorporated. The resulting problem is known as the (stationary) Stokes-Brinkman problem: 


T -I- div(— jxVw-f p/) =/ in D 

divM =0 in D 


( 2 ) 


We first introduce the finite element spaces, followed by the definition of the bilinear forms corresponding 
to the involved operators. 


2.1.1. H{div)-conforming Finite Elements for Stokes. 

Following m a DG formulation for the incompressible Navier Stokes equations which is locally conser¬ 
vative and energy-stable at the same time has to provide discrete solutions which are exactly divergence-free. 
This can be achieved with a suitable pair of finite element spaces. We consider the use of iJ(div)-conforming 
Finite Element spaces for the velocity u. H (div)-conformity requires that every discrete function Uh is in 

iF(div, D) = {u G [L 2 {D.)Y' : div u G L 2 {D)} 
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We use piecewise polynomials, so that on each element the functions are in i/(div, T). For global conformity, 
continuity of the normal component is necessary, resulting in 

Wh := {UT G n [^''(^)]^ luT-n]F = 0 W F G Th} d i?(div,0), Nw := dim(W,,), (3) 

TdTh 

with [•JiT’ the usual jump operator and the space of polynomials up to degree k. We refer to [H] for 
details on the construction of (div)-conforming finite elements such as the Brezzi-Douglas-Marini finite 
element. 

The appropriate Finite Element space for the pressure is the space of piecewise polynomials which are 
discontinuous and of one degree less: 

Qh-.= fVQ :=dim(Q^). (4) 

TdTh 

This velocity-pressure pair Wh/Qh fulfills 


AiY{Wh) = Qh- (5) 

The crucial point of this choice of the velocity-pressure pair is the property: 

If a velocity ut G Wh is weakly incompressible, it is also strongly incompressible: 


/ div{uT)q dx = 0 Vg/j G Qh ^ div(Mr) = 0 in fl. 
JQ 


( 6 ) 


The benefit of ([^ is twofold: First, it allows to show energy-stability for the incompressible Navier-Stokes 
equations, cf. section 2.5 Secondly, error estimates for the velocity error can be derived which are indepen¬ 


dent of the pressure field, we comment on this in remark below. 

As solutions of the incompressible Navier Stokes equations are x L 2 (^^)-regular and tangential 

continuity is not imposed as an essential condition on the finite element space the discrete formulation has 
to incorporate the tangential continuity weakly. We do this with a corresponding (Hybrid) DG formulation 
for the tangential components across element interfaces 


Remark 1 (Reduced spaces). Due to ([^ solutions of the discretized (Navier)-Stokes problem will he 
exactly divergence-free velocity fields. This a priori knowledge can he exploited. The basis for the space Wh 
can be constructed in such a way such that we can discard certain higher order basis functions (with non-zero 
divergence) that will have no contribution. A corresponding basis is introduced in 15If in the context of a 
set of higher order basis functions which fulfill an exact sequence property. The resulting space Wff^ has 
divfWJf^) = := {p\t = const : T G Th} such that also most of the degrees of freedom of the pressure 

space can be discarded. The reduction of basis functions for Wh and Qh is explained in more detail in ]2^ 
Chapter 2.2]. 


2.1.2. The HDG space for the velocity. 

To (weakly) enforce continuity we apply a discontinuous Galerkin (DG) formulation such as the Interior 
Penalty method m- However, to avoid the full coupling of degrees of freedom of neighboring elements, 
we introduce additional unknowns on the skeleton, the facet unknowns^ which represent an approximation 
of the tangential trace of the solution. The DG formulation is then replaced with a corresponding HDG 
formulation, s.t. degrees of freedom of neighboring elements couple only through the facet unknowns. We 
note that the resulting space is only “hybrid” in the sense of the tangential unknowns. Further, we do not 
consider a hybridization with respect to the pressure, cf. also remark 

As normal continuity is already implemented in Wh we only need a DG enforcement of continuity in the 
tangential direction and hence only introduce the facet unknown for the tangential direction of the trace: 

Fh ■= {UF G Yl [iP'=(F)]^ • n = 0 }, Nf:= diin{Fh). (7) 


5 




iJ^-conforming (CG) 



Figure 1: Tangential and normal continuity for different finite element spaces for the velocity. 


Functions in F/j have normal component zero. For the discretization of the velocity field we use the composite 
space 

Uh := Wh X Fh, Nu := dim{Uh). (8) 

Remark 2 (The role of the facet space). We note that the space Fh is only introduced to allow for a 
more efficient handling of linear systems. In section \4-t^ we consider a numerical example where a vector¬ 
valued Poisson problem is considered for the HDG space Uh in comparison to other (H(div)-conforming) DG 
spaces to illustrate the impact of Hybridization. In case that only explicit applications of discrete operators 
are used, the introduction of the space Fh entails no advantages. In this paper, facet variables appear only in 
the discretization of the viscous forces. Although the facet variable has no contribution in the discretization 
of inertia and pressure forces or the incompressibilty constraint we define the corresponding operations for 
u G Uh instead of ut G Wh to simplify the presentation. 


2.1.3. Viscous forces. 

For ease of presentation we consider a hybridized version of the Interior Penalty method for the dis¬ 
cretization of the viscous term. In remark we comment on alternatives. In the hybridized version of the 
Interior Penalty method, the usual jump across element interfaces of the Interior Penalty method is replaced 
with jumps between element interior and facet unknown (in tangential direction) |u*] = u^j, — up from both 
sides of a facet. The bilinear form corresponding to the HDG discretization of viscous forces is 


Ah{u, v) 



u = (ut,up),v=(vt,vf) G Uh 


(9) 


In this bilinear form the four terms have different functions. While the first two terms ensure consistency 
(in the sense of Galerkin orthogonality) the third and fourth term are tailored to ensure symmetry (adjoint 
consistency) and stability, respectively. Due to continuity of the solution (pj = 0) the latter terms also 
preserve consistency. For a more detailed introduction we refer to [221 section 2.3.1]. With respect to the 
discrete norm 


E 

TdTh 


||Vut||t + ^11 P*]llaT + ^ 


dup 

dn 


2 

dT 


which is an appropriately modified version of the discrete norm typically used in the analysis of Standard 
Interior Penalty methods, the bilinearform Ah is (for a sufficiently large a) consistent, bounded and coercive. 

The coupling through facet unknowns enforces the same kind of continuity as in the Interior Penalty 
DG formulation while preserving the following structure in the sparsity pattern: element unknowns are only 
coupled with unknowns associated with the same element or unknowns associated with aligned facets. 
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k 

4 

8 

16 

32 

quadrilateral 

0.305 

0.313 

0.315 

0.315 

triangle 

0.167 

0.190 

0.201 

0.205 


Table 1: LBB constant in dependence of k for a single element. 

Remark 3 (Interior Penalty and alternatives). A drawback of the (Hybrid) Interior Penalty method 
is the fact that the stabilization parameter a depends on the shape regularity of the mesh. Often a is chosen 
on the safe side, but as was pointed out in l,5gj / the condition number of arising linear systems increases 
with a. In the same paper a hybridized variant of the Bassi-Rebay stabilization method (cf. 0 0 [SBl) for 
a scalar Poisson equation has been proposed, see also remark^^ Such a variant is used in the numerical 
examples, but as it does not have any further consequences for the remainder of this work, we stick to the 
well-known (Hybrid) Interior Penalty method for ease of presentation. 

2 . 1 . 4 . Mass bilinear form. 

The HDG mass matrix is defined as 


M.)({u,v) := / UTVrdx, u = {ut,uf),v = {vt,vf) &Uh- 

Jn 

Note that we defined the mass matrix for u G Uh although uf has no contribution, cf. remarkj^ 

2.1.5. Pressure force and incompressibility constraint. 

For the pressure force and the incompressibility constraint we define the bilinearform 


( 10 ) 


for which the LBB-condition 


Vh{u,p) - j P divuT dx for u = {ut,uf) GUh, pG Qh 

T^h{u,p) 


T&n 


sup 

UGUh 




— ^LBB 


||p||l2, \/pGQh 


( 11 ) 


( 12 ) 


holds true for a independent of the mesh size h, cf. [HJ Proposition 2.3.5]. 

Remark 4 (Robustness in polynomial degree k). In numerical experiments we observed that the LBB 
constant c^^^ in (121 is also robust in k. To this end we computed the LBB constant of one element with 


Dirichlet boundary conditions and the condition f^pdx = 0 on the pressure. The results are shown in Table 
[7] /or varying the polynomial degree k. From the boundedness of the LBB constant on one element and the 
h-robustness in (12), the robustness in h and k on the global spaces follows. This will is in the forthcoming 


master’s thesis of Philip Lederer. 

2.1.6. Discretization of the Brinkman-Stokes problem 

With the introduced discretizations of the spatial operators the discrete Brinkman-Stokes problem can 
be written as: 

Find u G Ufi and p € Qh, s.t. 


r ^Mh{u,v)-{-Ah{u,v) + Dh(,v,p) = {f,v) 'ivGUh, 

T^h{u,q) = 0 V g G Qft,. 


(13) 


Due to coercivity of Ah (respectively r ^Mh+Ah), the LBB-condition of Vh and consistency and continuity 
of all bilinear forms Brezzi’s famous theorem (see m) can be applied to obtain optimal order a priori error 
estimates. For a discussion of the coupling structure of this discretization we refer to Remark below. 
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Remark 5 (Pressure-independence of velocity error). Classical error estimates for mixed problems 
result in error estimates which are formulated in the eompound norm of the velocity and pressure space. 
This has the disadvantage that the discretization error in the velocity depends on the approximation error 
in the pressure. As was pointed out in jhff, ideally this should not be the case. Due to the fact that discrete 


solutions to (131 are exactly divergence-free an error estimate for the velocity field which does not involve 
the pressure can he derived, cf. I2‘A Lemma 2.3.13]. 


2.2. Projected jumps: An enhancement of the HDG Stokes discretization 


The proposed HDG formulation for viscous forces in section 2.1.3 can also be derived as a hybrid mixed 
method with a modified flux. This is done in detail in [25] (see also [22| section 1.2.2] or [33]). In that 
setting the unknowns for the primal variable (wt) are approximated with the same polynomial degree k as 
the facet unknowns (mf)- Afterwards, in a postprocessing step approximations for the primal unknown, of 
one degree higher. A: + 1, are reconstructed in an element-by-element fashion. This approach ends up with 
a higher order approximation than the previously introduced HDG method considering the use of the same 
polynomial degree k for the facet unknowns. This sub-optimality can be overcome by means of a projection 
operator which leads to the projected jumps formulation. We first introduce the method and explain how 
the method can be implemented afterwards. 


2.2.1. Projected jumps: The method. 

The idea of the projected jumps formulation is to reduce the polynomial degree of the facet unknowns in 
Fh to k- 1, 

Fh^Fh:={uFG Y[[P'^-\F)f,UF-n = 0,\/FGFh}, (14) 

while keeping the polynomial degree k in Wh. In order to do this in a consistent fashion we have to modify 
the bilinearform Ah- First, we introduce the L 2 projection H for a fixed facet F € Fh'- 


H : [pfc(^)]d ^ [pfc-i(^)]d, f (^nf)vdx=[ fvdx ^ v&[V’^-^{F)Y (15) 

JF Jf 

Due to the fact that € V^~^{F) on affine linearly mapped elements, the test functions vt and vf in 
the first integral on the element boundary in Q can be replaced by their L 2 {F) projections Avt and Hwf- 
If we additionally use instead of |m*] to symmetrize and stabilize the formulation in (|^ we end up 

with a modified version of the Hybrid Interior Penalty method: 


Al{u,v) 



/ n|M*]n|?;‘] ds, 

JdT n 


u,v GWh X Fh 


(16) 


Note that the first two boundary integrals are just reformulated while only the last integral is really changed. 
This modification preserves the important properties of Ah- It is still consistent and bounded. Further, a 
deeper look into the coercivity proof (see [22]) reveals that coercivity can easily be shown in the modified 
(weaker) norm 


II?.* 



VutII 


diP“' 


lar' 


duT 


dn 


dT , 


(17) 


Notice that, as we do not modify Wh, the normal component of ut € Wh is still a polynomial of degree k. 

The idea of applying such a reduced stabilization is also discussed and analyzed in [snisD]. However, only 
for the two-dimensional case the realization of the projection H is discussed (by means of Gauss quadrature, 
cf. [391 section 3.4]). In the next section a simple way to implement this projection operator is presented. 









Remark 6 (Interplay with operator-splitting). If a Hybrid DG formulation for a problem involving 
diffusion (viscosity) and convection is applied, such a reduction of the facet unknowns is only appropriate if 
diffusion is dominating. One premise of the operator-splitting in this paper is that convection is not involved 
in implicit solution steps. Hence, the facet variable will never be involved in the discretization of convection, 
s.t. even if the physical problem is convection dominated, the projected jumps modification can be applied 
without loss of accuracy. 


Remark 7 (Comparison to a fully hybridized formulation). We note that the formulation in (131 
(and the corresponding reduced formulation with is not a “hybridized” formulation in the usual 

sense, see e.g. f25^ . Only for the tangential component of the velocity a new unknown is introduced, such 
that we only have a “tangentially hybridized” velocity discretization. The pressure is not hybridized. We 
comment on the resulting coupling of this HDG method and compare it to a HDG formulation where only 
Hybrid unknowns on the facet exist, for the velocity and the pressure. The velocity unknowns inUh = Wh x 
(or Wh X Fh in the reduced case) can be reduced to unknowns on the skeleton by static condensation. The 
remaining facet unknowns are the usual unknowns for the normal component of the H{div)-conforming 
finite element space and the (tangential) unknowns in Fh (or Fh). Due to the pressure unknowns, static 
condensation of the formulation in (13) does not yield a global system which only involves facet unknowns, 
but also element unknowns for the pressure. However, except for the constant pressure all pressure unknowns 
can be eliminated by static condensation. We summarize: The globally reduced system only involves velocity 
unknowns on the facets and a constant pressure on each element. Alternatively one could formulate a 
related HDG formulation which only involves facet unknowns for the velocity and the pressure. In such a 
formulation the pressure plays the role of the lagrange multiplier for the normal continuity and only for 
the (weak) tangential continuity velocity unknowns have to be added on the facet. To obtain an accurate 
(and exactly divergence free) order k approximation of the velocity and an order k — 1 approximation of the 
pressure, the facet unknowns would be an order k approximation of the pressure and an order k (or order k—1 
if projected jumps or corresponding postprocessing techniques are applied) approximation of the tangential 
velocities. Overall the number of unknowns would be smaller compared to the formulation proposed before 
by one (constant) pressure unknown per element. However, the price for this is the fact, that the Lagrange 
multiplier space (the pressure unknowns on the facets) is much larger and hence the structure of the arising 
linear system is different. It very much depends on the linear solver strategies which of the two Hybrid DG 
formulations gives the better perfomance. 


2.2.2. Projected jumps: Realization. 

The following way to implement projected jumps needed for Aj, in (16) relies on an L 2 -orthogonal basis 
for the facet functions in Ff,. To obtain an L 2 “Orthogonal basis, we take a local coordinate system on 
each facet (spanned by d — 1 aligned edges) and an L^-orthogonal Dubiner basis (see [3S]) for each vector 
component. We consider the three dimensional case, and denote the vectors of the local coordinate system 
by ei and 62 - Translation to the two dimensional case is then obvious. 

On each facet F G Fh, the space of (vector-valued) polynomials up to degree k, span(ei, 62 ) • V^{F), can 
be split into the orthogonal subspaces 


Vfe_i := span{ei,e 2 } • T’'" ^(F) and := span{ei, 62 } • [P''(F) n (T’'" ^(^))'^]- (18) 

On each facet F the facet unknown up can thus be written as up = up + i^ with unique up G Vfe_i, A G Vh-i- 
We now replace the highest order function A with two copies At and At' , each of these functions is associated 
with one of the two neighbouring elements T and T', the functions At are only defined element-local. At 
can thus be eliminated after the computation of the element matrix and finally, only up G Wh and up £ Fh 
appear in the global system. 

At is only responsible for implicitly realising the projection operator. To see this we now consider one 
facet F £ Fh and one of the neighboring elements T GTh. We use the decomposition corresponding to ( [I^ 
for trial and test functions 

Up = Up At and vp = vp pp, (19) 
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with uf,vf G Vfc_i and Xt^Ht G where are only supported on T. In (13) we choose the test 

function v such that vt = vf = 0, G V^_i on T and /ir' = 0 on every other element T'. This yields 


[ ^ ds + [ v—{v!g — uf — ds = [ — At)mt ds = {), V /iy € V, 

Jf on Jf h Jf h 


t-i- 


( 20 ) 


Hence, there holds Xt = {I — Il)u^ and we have 


[u*] = Uf — Uf — Xt = n(M^) — uf = n(Mp — uf) = n|r 


( 21 ) 


We conclude that it is sufficient to consider the HDG bilinear form Ah as before, where the local element 
matrices are computed according to Fh- For each element matrix the degrees of freedom corresponding to 
"^k-i then eliminated forming a corresponding Schur complement. This yiels a final stiffness matrix 
which is only set up with respect to space unknowns of Wh y. Fh- 


2.3. DG formulation for the convection 

For the discretization of the convection part we consider a standard DG finite element space: 


Vh:={u-.UG [V'^{T)Y vt G r4, Nv := dim(ld,). 


We define the mass bilinear form 



w z dx, 


w,z G Vh- 


( 22 ) 


(23) 


Using an L^-orthogonal basis on each element renders the associated mass matrix : RW RW 
diagonal. A stable spatial discretization of ([^ with respect to the convection part is achieved using a 
Standard Upwind DG trilinearform. We assume that the convection velocity ut is exactly divergence-free. 

Ch{uT',w, z) / w^UT'.yz dx + 

-r Jt 


UT-n w z ds, Ut G Wh, w,z G Vh- 


IdT 


(24) 


where w denotes the upwind value w = lirngx^o vj(x — eut^x)). Standard Upwind DG formulations are stable 
in the sense that with := {x G ut{x) • u < 0} there holds 

- j \ut ■ n\ ds + Ch{uT', w,w) > 0, V w G 14, u G Uh, div(Mr) = 0. (25) 

2 Jaair, 


2-4- Transfer operations and embeddings 

The discrete convection and Stokes operators have been defined on different spaces. In order to combine 
both we introduce transfer operations between the (finite element) spaces to make both discretizations 
compatible. We restrict ourselves to two types of transfer operations which are based on embeddings: 

I: With Wh C 14 we have a canonical embedding of Uh in 14 with the embedding operator X : (ut, uf) G 
Uh ^ Ut G Vh- Note that the corresponding operation in terms of coefficient vectors, denoted by 
I : R^^ —>• R-^v ^ jg identity. 

I^: The embedding operator X implies the canonical embedding X' : Vf ^ U'^, which maps functionals 
I' : / G 14 —)■ [{ut, Uf) G Uh —>■ /(mt)] G Uf. The corresponding matrix representation is : R'^'^ —)■ 


To realize the operator I we consider the equivalent problem for u G Uh- 




In terms of coefficient vectors this reads as 


M^(Iu) = Vue 

with the mixed mass matrix 


I = 


(27) 




(fiY^pY dx, i = = Nw and = 0, j > Nw 


uy 


Note that is diagonal (for affine linear transformations) such that (M^)“^ can be evaluated very 
efficiently. The overall cost of the transfer operator I is essentially that of one sparse matrix multiplication. 

Let x —>• denote the discrete convection operator corresponding to the trilinearform 

Ch in ( |2^ . 

With the operator I we can formulate applications of the convection and the mass operations C^(u), : 

—)■ with respect to functions in the HDG space Uh and denote the corresponding operators by 


C^(u) := l'^C''(u)I, := I^M^I. (28) 

Remark 8 (Restriction on time integration scheme). Note that the restriction to these two transfer 
operations implies that we do not allow to apply any part of the Stokes operator to a function in Vh and that 
no functional on Uh can be used in solution steps involving the convection operator. This is a restriction on 
the time integration scheme. 

Remark 9 (Curved elements). If curved elements are considered we no longer have Wh C Vh due to the 
Piola transform usually applied to construct H{div)-conforming finite elements. In this case I : Uh ^ Vh 
is not an embedding, but the projection. Nevertheless, is block diagonal with blocks which are not 
diagonal matrices only on curved elements. Hence, applications o/(M^)“^ are still cheap. 


2.5. The semidiscrete formulation 

With the definitions of the bilinear forms in the previous sections we arrive at the following spatially 
discrete DAE problem: Find u{t) G Uh and p{t) G Qh, such that 

f MY,{§iU,v)+Ah(u,v)+Ch(uT-,u,v)+Vh{v,p)= {f,v) \/v e Uh, tG[0,r], 

Vh{u,q) = 0 yq€Qh, tG[0,T], (29) 

[ Mh{u,v)=Mh{uo,v) Wv^Uh, t = 0. 


Here, we implicitly used the embedding T to define Ch{uT', ■,■) on Uh. Due to (25) we have with 


i^U,u)L2 = = hh^^Wuh^ = {f,u) - Ah{u,u) -Ch{uT-,U,u) -Vh{u,p) < ||/||i2||u||i2 (30) 


>0 


>0 


=0 


the stability of the kinetic energy: 


<||/(t)||L. 


(31) 


In the next section we discuss operator splitting time integration methods to solve (29) efficiently. 


3. Operator-splitting time integration 

In this section we are faced with the problem of solving the semi-discrete Navier-Stokes problem with a 
proper time integration scheme. For ease of presentation we neglect external forces (f = 0) in the following 
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and consider the problem in terms of discrete operators corresponding to the bilinear (trilinear) forms 
introduced in the previous section: Find u{t) G and p(t) G such that 


+ Au +C^(u) u+AtDp = f in [0,T], 
D^u = 0 in [0,T], 
u(t = 0) = Uq. 


(32) 


In the time integration scheme we want to explicitly exploit the properties of the spatial discretization, 
i.e. the convection operator C should only be involved explicitly in terms of operator evaluation. Due 
to the DAE-structure we require time integration schemes which are stijfly accurate. Hence, the solution 
of a Stokes-Brinkman problem should conclude every time step so that the incompressibility constraint is 
ensured. For operator splittings of this kind different approaches exist. We briefly discuss three approaches. 


In section 3.1 we discuss additive decomposition methods, like the famous class of IMplicit EXplicit (IMEX) 


schemes. Product decomposition methods, sometimes also called exponential factor splittings, are discussed 

We discuss the 


3.2 in the framework of operator-integration-factor splittings introduced in 


in section 

advantages and disadvantages of both approaches and motivate the consideration of a different approach. 
An operator-splitting modification of the famous fractional step method, cf. |56j . which eliminates the 
most important disadvantages of the additive and multiplicative decomposition methods is then introduced 
and discussed in section |3.3[ We want to stress that the considered operator splitting approaches are of 
convection-diffusion type and should not be confused with projection methods like the Chorin splitting EZ]. 

3.1. Additive decomposition using IMEX schemes 

Additive decomposition methods distinguish spatial operators that are treated implicitly and those that 
are treated explicitly. The decomposition is additive in the sense that every solution (sub-)step involves 
both operators. This is used by IMplicit EXplicit (IMEX) schemes (see [HI EHl ES] )■ The simplest of these 
schemes is the semi-implicit Euler method for which one time step of size At reads as: 


(M'^ -I- AtA)u"+i4- AtDp"+i = - AtC'^(u")u” 

D^u^+i = 0 


(33) 


Convection only appears on the r.h.s. so that only operator evaluations occur for the convection operator, 
while the remainder appears fully implicit. This scheme is obviously only first order accurate and only con¬ 
ditionally stable. For higher order methods of this decomposition type (e.g. using Multistep or partitioned 
Runge-Kutta schemes) we refer to [HI 03 05] . 

The major drawback of this type of operator splitting methods is the fact that the time step size of the 
explicit and the implicit part of the decomposition have to coincide. Thereby the stability restriction caused 
by the explicit treatment of the convection part dictates not only the number of explicit evaluations but also 
- which is typically more expensive - the number of solution steps for the implicit part. The decomposition 
methods considered in the sections |3]2] and |3l3] overcome this issue. 


3.2. Product decomposition with operator-integration-factor splitting 

The disadvantage of IMEX methods can be avoided with product decomposition methods, where se¬ 
quences of separated problems are solved successively. The separated problems then only involve the Stokes 
or the convection operator at the same time. The major benefit of this is the fact, that the numerical solu¬ 
tion of the sub-problems (e.g. the size of the time steps) can be chosen completely different. The derivation 
of the so called operator-integration-factor splitting approach has been derived in |48l Section 2.1]. We use 
this idea to obtain an operator decoupling for the DAE which allows to formally rewrite (32) as 


dt 


(Q*^**u)-hg‘^**M-i(Au + Dp) =0 


= U m 


D^u =0 


[0,T], 

[0,T], 


(34) 


for some arbitrary t* 
section 13.2.21 


with the integration factor , the propagation operator, specified later in 
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Be aware that M ^ is to be understood only formally for now. We are able to apply any suitable (i.e. 


implicit and stiffly accurate) time integration method for (34). We do this for a first order method here to 


explain the procedure and refer to the literature [48] for higher order variants. 
3.2.1. Implicit Euler. 


We employ the implicit Euler method on (34) and arrive at 




After setting t* = and multiplication with M this simplifies to 

(M +AtA)u"+i + AtDp"+i) 

D^u"+i = 0, 


-Dp"+i) =0, 

D^u"+i = 0. 


lu" 


(35) 


(36) 


which is the solution of a Stokes-Brinkman problem as in (33), but with a different right hand side. 


3.2.2. The propagation operator. 

The propagation operator Q* is defined as Q* 


= v(t 2 ) where v(s) solves 


9v 

Ws 


+ (M^) (s)v(s) = 0, VsG(tyr], v(ti) =-w. 


(37) 


Here C^(s) := C^(u(s)) with u(s) an extrapolation of divergence-free solutions of previous time steps. 
Note, that the extrapolation of convection velocities renders the problem (37) linear hyperbolic and ensures 
stability in the sense of (p^. After replacing Q by a numerical time integrator Qa* for (37), the time 


integration method is completely specified. The order of accuracy of the extrapolation u and the time 


integrator Qa* should coincide with the order of the time integration method applied on (34). 


3.2.3. Properties of product decompositions. 

At this point the advantage of product decomposition methods over IMEX schemes is evident: As the 


time integration of (37) is independent of the one for (34) the (not necessarily) explicit time integrator can 


deal with stability restrictions (typically using multiple time steps) without influencing the time step for 


(341. This separation of the two problems also introduces the biggest disadvantage of product decomposition 
methods: an additional consistency error. Even if the DAE (32) has a stable stationary solution, product 
decomposition methods may not reach it. This is not the case for monolithic or additive decomposition 
approaches. Nevertheless, this splitting error is controlled by the time discretization error. 


Remark 10 ((Marchuk-)Yanenko splitting). If the implicit Euler discretization as in (36) is combined 


with an Euler method for (37), the famous Yanenko splitting method is recovered. 


3.3. A modified fractional-step-9-scheme 

In this section we introduce an approach to circumvent severe CFL-restrictions and splitting errors at 
the same time. We no longer ask for sub-problems that only involve the Stokes or the convection operator 


(as in section 3.2), but ask for sub-problems which only involve one of both implicitly. The method is based 
on [33] where the well-known fractional-step-0-scheme, cf. (SHI Chapter II, section 10], is modified. The 
resulting method is an additive decomposition method without the time step restrictions of IMEX schemes. 

We start with formally writing down an operator-splitting version of the fractional-step-0-scheme. The 
scheme is divided into three steps, the first and the last step treat the Stokes part implicitly and the 
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convection part explicitly as in (331 while the second step treats convection implicitly and viscosity forces 
explicitly. 


Step 
Step 2 
Step 3 

The time stages are labeled by the superscripts no, ni, n 2 , n^, respectively, where no denotes initial data 
and ns the final time stage. The time steps size for the first and the last step is 9At and 9*At in the 
middle step with 9* = 1 — 29. Here 9 = 1 — l/-\/2. We note that specifications for M and C with respect 
to the considered spaces and the convection velocity for C are still missing at this points. The scheme is 
second order accurate. In contrast to the unsplit fractional-step-0-scheme the stability analysis of this time 
integration scheme is an open problem. In our experience, however, the stability restrictions are much less 
restrictive than those of a comparable IMEX schemes. 


, (M + 6»AtA)u”i+6»AtDp”i =Mu"'>-6»AtCu"o 

^ = 0 


r=):{ (M + rAtC)u” 


= Mu”i - 6»*At(Au”i + Dp”i) 


„3. , (M + 6»AtA)u”3 + 9At D p”" = - 9AtCu^^ 

* ■ 1 = 0 


(38a) 

(38b) 

(38c) 


3.3.1. Sub-steps in different spaces. 

Initially, we stated that we want to avoid solving linear systems involving convection, for efficiency 


reasons. This seems to be contradictory to what is formulated in (38b). Nevertheless, as only convection is 


involved implicitly an efficient numerical solution is still possible. We apply a simple iterative scheme which 
only involves explicit operator evaluations of the convection to do so. We explain this in section |3.3.2[ To 
do this efficiently, we want Step 2 to be formulated in the space 14, while Step 1 and 3 are to be formulated 
in the space {Uh, Qh)- This poses problems the solution of which we discuss in this section. 

In Step 1 and Step 3 the adjustments are obvious: Step I only depends on initial data in Uh- The 
initial data for Step 3 is in 14 but appears only in terms of functionals (Mv and Cv) such that the transfer 
operations are clear. Step 2 is more involved, cf. remark]^ Functionals in (such as Au) are in general 
not functionals in Vff We use the first equation in ( |38a| ) to formally define a different representation of the 
functionals required in Step 2: 


9At = -6»At(Au”i + Dp”4 = M(u"i - u”“) + 9AtCu^° 

Now we replace M and C with operations suitable for a setting in 14: 


9 g"i := — - u”“) + 6lAtC^(u"“)Iu"« 

9At ^ ' ’ 

We arrive at the modified fractional-step-0-scheme: 

/ (M'^+ 6»AtA)u’"i+ 9AfD p’"i = - 9At l'^C'^(u”«) I u"o 

btep 1 : < ^ 0 


Step 2 : { (M^'t 9*AtC^{t^^)^ 


= - 6»*Atg" 


Step 3 : 


(M'^+ 6»AtA)u”3+ eAtT> p"3 = 9At 

D'^u”3 = 0 


(39) 

(40) 

(41a) 

(41b) 

(41c) 


where we replaced the generic operators M, C with suitable ones. We recall the definition of the extrapolated 
convection operator C^(t"2) := C'^(u(t”^)) the convection velocity of which is (linearly) extrapolated from 
exactly divergence-free velocities. 

3.3.2. Iterative solution of the implicit convection problem. 


In (41b) we need to solve a problem of the form 

v-kT*(M'^)-iC'^v = g* 


14 


(42) 








for given r*, g* and constant convection C^. We do this by means of a pseudo time-stepping method, i.e. 
we formulate (421 as the stationary solution to 


^v(s) +v(s) + T*C^v(s) = g*, se[0, oo). (43) 

Note that a stationary solution exists as Id + r*C'^ only has eigenvalues A, with Re{X) > 1. This stationary 
solution is approximated with a few explicit Euler time step with an artificial time step size As. 

Vi+i =v, + As(g* -\'i + (44) 


This procedure can be interpreted as a Richardson iteration applied to (42). With a time step size which is 
tailored for stability (as in the numerical solution to ( [37| )) we iterate (44) until the initial residual is reduced 
by a prescribed factor. To solve the convection step, Step 2, we only require operator evaluations for the 
convection. The time step size As in this iteration is decoupled from the time step size At used of the 
overall scheme. 


4. Numerical examples 


In this section we consider different test problems which essentially purpose three different goals: 

1. The validation of convergence properties of the HDG method with and without the projected jumps. 

2. The investigation and quantification of the dependency of the sparsity pattern of arising linear systems 
on the choice of the discrete velocity spaces. 

3. The evaluation of the performance of the space and time discretization on benchmark problem. 


The examples in section |4.1| and section |4.2| aim at goal 1 and are two-dimensional test cases with known 
exact solutions for a stationary Stokes and Navier-Stokes problem, respectively. These cases allow for a 
thorough investigation of the convergence history of the proposed spatial discretizations in the usual norms. 
The test case in section |4.3| aims at goal 2 and is concerned with the impact of the choice of discretization 
spaces on the complexity of linear systems. For this purpose, a three-dimensional vector-valued Poisson 
problem is considered using different i?(div)-conforming DG and HDG spaces. The impact of hybridization 
and the improvement due to the projected jumps modification is compared to other DG methods. Finally, 
we approach goal 3 by considering two challenging transient benchmark problems in sections |4.4| and |4.5[ 
The benchmark problems have been defined within the DFG Priority Research Programm ’Flow Simulation 
on High Performance Computers’ and are formulated in |58j . In section 4.4 a two-dimensional benchmark 


problem is considered to demonstrate the accuracy of our method for a demanding test case and to compare 
the discussed time integration methods. Finally, in section |4.5| we discuss a three-dimensional, and hence 
computationally demanding, benchmark problem from [58] . We compare accuracy and run-time performance 
to the data of the studies in 1201 • 

The methods discussed in this paper have been implemented in the add-on package ngsf low m for the 
high order finite element library NGSolve [62]. The computations in this section have also been carried out 
with this software. Throughout this section we only consider direct solvers and comment on linear solvers 
below, in remark 1121 


For the computations of the benchmark problem in sections 4.4 and 4.5 we used the reduced H{div)- 

red 


conforming space cf. remarkand the projected jumps modification, cf. section 2.2 


4-1- Stokes flow around obstacle 

We consider the Stokes problem in the domain D = [—2, 2]^ \ il~ with the circular obstacle := 
{ll^^ll < 1} and viscosity v = 1. On the whole boundary except for the outflow boundary {x = 2} we 
prescribe Dirichlet boundary conditions, on {x = 2} we prescribe Neumann-kind boundary conditions such 
that the solution to the problem is given by u = (Sy^/, —clj,^') with the potential 'I' = y{l — ^2^y2 )■ Note 
that 'I'lan- is constant so that u ■ n = V'k x n = 0. On we have u ■ n = 0, but u ^ 0, i.e. we have a 
slip on the obstacle. Further we have that p = 0 in D. 
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Figure 2: Convergence of divergence-conforming HDG method with (solid) and without (dotted) projected jumps in different 
norms for the example in section [4.1| The dotted lines in gray correspond to the HDG method without projected jumps, the 
solid lines in color correspond to the HDG method with projected jumps. Except for the error in the pressure the results are 
hardly distinguishable. 


On an initially coarse unstructured grid with only 72 triangles we obtain the results displayed in Figure 
[Rafter succesive uniform mesh refinements. We compare the exact solution u, p and the computed solution 
Uh, Ph and investigate the convergence of the error in the norms \\u — Uh\\L^(Q), ||V(w — and 

\\p ~ Ph\\L^{n)- We consider the discretization with and without the projected jumps modification. In all 
norms we observe optimal order convergence, i.e. \\u — Uh\\L^(Q) = ||V(u — = 0 {h^) and 

lb ~ Ph\\L^(ci) = 0{h^) where k is the order of the velocity field inside the elements. Note that after static 
condensation the remaining degrees of freedoms are the unknowns corresponding to order k polynomials 
for the tangential and order k polynomials for the normal component on the facets and one constant per 
element for the pressure. For the projected jumps formulation the tangential unknowns are reduced by one 
order. The difference between the error of both formulations — with and without projected jumps — is only 
marginal in the velocity. In the pressure field we observe a difference between the methods, however both 
converge with optimal order and the difference decreases for increasing polynomial degree k. For the impact 


of the projected jumps modification concerning the computational effort, we refer to section 4.3 where this 
aspect is discussed in more detail. 


4-2. Kovasznay flow 

As a test case for the stationary Navier-Stokes equations, we consider the famous example by [63]. On 
the boundary of the domain O = [—^ 4]x [0, 2] we again prescribe inhomogeneous Dirichlet data, so that 
the exact solution to the Navier-Stokes equations (with z/ = I) is 


L{x,y) = 


1 — cos(27rj/) 
sin(27r?/) 


p{x,y) = --e 


2\x 


' p with p G 


Here, A = 




- 2 _)_g 4^'2 s-iid p SO that f^pdx = 0. 

To approximate the stationary solution we initially solve the Stokes problem corresponding to the bound- 


+\Ay 


ary data and use the simple IMEX scheme of first order, cf. (33), to progress to t = 1000 with 50 time steps 


of size At = 20. This choice of time discretization parameters gives stable solutions for all considered spatial 
discretizations but at the same time provides time discretization errors which are neglegible compared to 
the spatial errors. On an initally coarse unstructured grid with only 18 triangles we obtain the results 
displayed in Figurej^after succesive uniform mesh refinements. We observe that the errors converge optimal 
in the considered norms for the pressure and the velocity. Further, we observe that the difference between 
the results obtained with and without the projected jumps formulation, e.g. with and without a reduction 
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Figure 3: Convergence of divergence-conforming HDG method with (solid) and without (dotted) projected jumps in different 
norms for the example in section f4.2| The results are almost identical. 
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refinement level 


of the degrees of freedoms at the facets, is again only marginal. In fact one only observes a (very small) 
difference in the case k = 1. 


4-3. Linear systems - A comparison between DG and HDG methods 
We consider the comparably simple vector-valued Poisson problem 


— Am = / in fl, u = 0 on dLl, 


(45) 


discretized with four different methods. The three-dimensional domain LI is the same as the one in section 14^ 
with an unstructured tetrahedral mesh consisting of 3487 elements. The problem (451 leads to discretizations 
with symmetric positive definite matrices A. Only the lower triangular part of the sparse matrix has to be 
stored and a sparse Cholesky factorization algorithm can be used. In this comparison, we used the sparse 
direct solver PARDISO, cf. [64l[65]. We are only concerned with the sparsity pattern of different methods 
before and after static condensation and do not compare accuracy or conditioning of the discretizations. 
Moreover, for the sake of simplicity we do not apply the reduction of the space Wh as discussed in remark 
The following quantities of interest for the linear systems arising from discretizations of (451 are displayed in 
^dof[K] : number of unknowns (in thousands) 

^cdof[K] : number of unknowns after static condensation (in thousands) 


Table 


^nzeA[K] : number of non-zero entries in the system matrix A (lower triangle) (in thousands) 
^nzeL[K] : number of non-zero entries in the Cholesky factor L (in thousands) 


Four different methods are considered on Wh or Uh with varying polynomial degree between fc = 1 and 
fc = 6: 

1. HDG: The HDG method proposed in section without projected jumps. 

2. PHDG: The HDG method with projected jumps, cf. section 

3. Std.DG.: A standard DG method using the space Wh where the basis is constructed such that all 
degrees of freedom from one element couple with all degrees of freedom from adjacent elements. 

4. N.DG: A nodal DG method using the space Wh where basis functions are assumed to be constructed 
such that degrees of freedom associated to one element couple only with degrees of freedom from 
adjacent elements which have support on the shared facet. At the same time we assume that basis 
functions are constructed such that the number of basis functions with support on a facets is minimized, 
cf. m- In terms of the sparsity pattern this nodal DG method represents the best case for a DG 
method without hybridization. 
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Std.DG 

N.DG 

HDG 

PHDG 

Std.DG 

N.DG 

HDG 

PHDG 



k 

^ 1 



k 

^ 2 


#dof[K] 

23 

23 

69 

38 

67 

67 

158 

112 

#cdof[K] 

23 

23 

69 

38 

67 

67 

137 

91 

^nzeA[K] 

732 

676 

2 037 

637 

5 073 

4 177 

8 103 

3 628 

^nzeL[K] 

10 113 

10 208 

17 768 

5 569 

64 463 

69 831 

70 138 

31 415 



k 

^ 3 



k 

^ 4 


#dof[K] 

146 

146 

298 

237 

271 

271 

500 

423 

#cdof[K] 

146 

146 

229 

168 

271 

261 

343 

267 

^nzeA[K] 

21 686 

16 051 

22 443 

12 124 

69 525 

46 912 

50 412 

30 588 

^nzeL[K] 

261 977 

260 416 

194 524 

104 496 

814 168 

731 253 

435 129 

264 183 



k 

^ 5 



k 

^ 6 


#dof[K] 

453 

453 

773 

681 

702 

702 

1 128 

1 022 

#cdof[K] 

453 

411 

480 

389 

702 

597 

640 

533 

^nzeA[K] 

184 035 

101 473 

98 696 

64 816 

424 764 

200 813 

175 321 

121 944 

^nzeL[K] 

2 099 690 

1 741 072 

847 913 

557 798 

4 752 072 

3 519 061 

1 502 558 

1 045 444 


Table 2: Comparison of different DG methods for the vector-valued reaction Poisson problem ||45|l 


In Table we observe that for small polynomial degree k the amount of additional unknowns required for 
the HDG formulation is quite large as are the nonzero entries in the system matrix. Nevertheless except 
for fc = 1, the number of nonzero entries in the Cholesky factor are comparable (fc = 2) to the Standard 
DG methods or less (fc > 2). For high order, i.e. fc > 4 the HDG method performs signihcantly better as it 
has less nonzero entries in A and L. The HDG space with the projected jumps modification improves the 
situation dramatically. It essentially compensates the overhead of the HDG method for small k. But even 
for fc = 6 the effect is still significant. We note, that the difference between the first three methods increases 
for increasing polynomial degree fc whereas the difference between the last two method, PHDG and HDG, 
decreases. In all cases the HDG method with projected jumps outperforms all alternatives. 

j.j. A two-dimensional benchmark problem 

In this section we consider the benchmark problem denoted as “2D-2Z” in [58] where a laminar flow 
around a circle-shaped obstacle is considered. The Reynolds number is moderately high {Re = 100) and 
results in a periodic vortex street behind the obstacle. We briefly introduce the problem (for more details 
we refer to [58] 1 and the numerical setup to investigate spatial and temporal discretization errors. Finally, 
we discuss the obtained results. 


4 . 4 . 1 . Geometrical setup and boundary conditions. 

The domain is a rectangular channel without an almost vertically centered circular obstacle, cf. Figure 

El 

n := [0, 2.2] x[0,0.41] \{||x- (0.2,0.2)112 < 0.05}. (46) 

The boundary is decomposed into Tin '■= {x = 0}, the inflow boundary. Tout ■= {x = 2.2}, the outflow 
boundary and Fpi/ := dTt\ (ri„UFo„t), the wall boundary. On Tout we prescribe natural boundary conditions 
{—vVu + pi) • n = 0, on Fw homogeneous Dirichlet boundary conditions for the velocity and on F^ the 
inflow Dirichlet boundary conditions 

m(0, j/, t) =ud = (3/2 • u) • 4 • y{dy - y)/dl ■ (1, 0,0). 

Here, u = 1 and the viscosity is fixed to = 10“^ which results in a Reynolds number Re = 100. 


Drag and Lift. 

The quantities of interest in this example are the (maximal and minimal) drag and lift forces cd, cl that 
act on the disc. These are defined as 


1 


CD '■ = 


du 

- 

on 


ds, 


1 


Cl ■= 

u^r 


du 
- 

on 


ds. 
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Here ex,ey denote the unit vectors in x and y direction, r = 0.05 is the radius of the obstacle, u is the 
average inflow velocity (u = 1) and To denotes the surface of the obstacle. 

4.4-3. Numerical setup. 

We use an unstructured triangular grid with an additional layer of quadrilaterals around the disk which 
is anisotropically refined towards the disk once. In Figure]^ the geometry, the mesh and a typical solution 
is depicted. 


2 


0 

Figure 4: Sketch of the mesh and the solution (color coding corresponding to velocity magnitude ||u|| 2 ) to the problem considered 
in section [ 4 .4| at a fixed time t (left) and zoom-in on the boundary layer mesh (right). 

In order to be able to neglect time discretization errors, when investigating the spatial accuracy, we 
consider the use of a well-known stiffly accurate second order Runge-Kutta-IMEX scheme, taken from [321 
section 2.6], with an extremely small time step size 3.125 • 10“® which means that roughly 10000 time steps 
are used to resolve one full period. The duration of one full cycle is roughly l/3s. We also fix the mesh and 
consider a pure p-refinement, i.e. variations of the polynomial degree k. 

For the investigations of the temporal discretization error, we consider a fixed polynomial degree fc = 6, 
with the same mesh. For the considered example the stability restriction of the second order IMEX scheme is 
severe. A time step of below 10“^ has to be considered. The intention of the discussion of operator-splitting 
methods in section has been to present strategies to circumvent or relax these severe conditions. We 
compare the performance of a product decomposition method, a second order operator-integration-factor 
splitting version of the BDF2 method, and the modified fractional-step-0-scheme discussed in section |3.3[ 
Note that these methods allows to consider much larger time steps than the IMEX scheme. As references 
we give the values from the literature [521 and from the IMEX scheme with At = 10“^ (stability limit) and 
the reference solution with At = 3.125 • 10“®. 

4.4-4- Numerical results: spatial discretization. 

In Table [^ the quantities of interest are shown for varying polynomial degree k. As a reference we also 
show the result obtained by FEATFLOW m with a discretization using quadrilateral meshes and continuous 
second order finite elements for the velocity with a discontinuous piecewise linear pressure {Q 2 /Pi'‘^'^). These 
results have been made accessible on [66]. We observe a rapid convergence for the p-refinement, i.e. the 



#dof 

max CD 

min CD 

max Cl 

min Cl 

k = 1 

2 211 

2.52594 

2.47871 

0.65728 

-0.81672 

fe = 2 

4 148 

3.22841 

3.16260 

1.00571 

-1.03894 

fc = 3 

6 558 

3.23184 

3.16842 

0.98822 

-1.02427 

k = 4 

9 441 

3.22714 

3.16401 

0.98431 

-1.01906 

k = 5 

12 797 

3.22759 

3.16432 

0.98578 

-1.02053 

II 

16 626 

3.22757 

3.16430 

0.98580 

-1.02053 

ref. [66] 

167 232 

667 264 

3.22662 

3.22711 

3.16351 

3.16426 

0.98620 

0.98658 

-1.02093 

-1.02129 




Table 3: Accuracy of the spatial discretization: results for different polynomial degrees. 


increase of the polynomial degree. Compared to the results from the literature [66] the same order of 
accuracy is achieved with a lot less degrees of freedoms. 
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4.4-5. Numerical results: temporal diseretization. 

The results for the temporal discretization are shown in Table First of all, we observe that the second 
order IMEX scheme is already very accurate at its stability limit. But the method does not allow to choose 
larger time steps. This is in contrast to the alternative methods discussed here. For these methods we 
can consider much large time steps and observe a second order convergence. In this example the product 
decomposition method is more accurate than the modified fractional step method by one (time) level. We 
note that one major concern with this method is however, that splitting errors appear also if a stationary to 
the flow problem exists. This is not the case for additive decomposition methods or the proposed modified 
fractional step method. 



1/At 

max CD 

min CD 

max Cl 

min Cl 


<1 000 


unstable 


2nd order IMEX 

1 000 

3.22754 

3.16437 

0.98486 

-1.01957 


32 000 

3.22757 

3.16431 

0.98580 

-1.02053 

operator-integration- 
factor splitting 
(BDF2) 

125 

3.38456 

3.28279 

1.16052 

-1.19939 

250 

3.25564 

3.18643 

1.01725 

-1.05276 

500 

3.23239 

3.16836 

0.98864 

-1.02378 

1 000 

3.22819 

3.16394 

0.98320 

-1.02013 

modified fractional- 

125 

3.38272 

3.29673 

1.07667 

-1.12227 

step-0-scheme 

250 

3.28713 

3.21483 

1.00937 

-1.04840 

500 

3.25036 

3.18127 

0.99172 

-1.02807 


1 000 

3.23656 

3.17094 

0.98721 

-1.02227 


Table 4: Accuracy of operator-splitting time integration methods 


4 . 5 . A three-dimensional benchmark problem 

Finally, we consider a three dimensional benchmark problem, the problem “3D-3Z” in [58]. In contrast 
to the problem in section [4.4| the inflow velocity is varied over time and the observed time interval is fixed. 
The maximal Reynolds number in this configuration is also Re = 100 as in the previous section. The focus 
in this section is on the study of performance in the sense of computational effort over accuracy. 


4-5.1. Geometrical setup and boundary conditions. 

The geometrical setup is a generalization of the problem in the previous section. A cylindrical-shaped 
obstacle is places in a cuboid-shaped channel slightly above the vertical center: 

b!:= [0,2.5]x [0,0.41] X [0,0.41] \{||(a;i, era) -(0.5,0.2)112 <0.05}. (47) 

Boundary conditions are chosen as in the previous example except for a change to unsteady inflow boundary 
conditions: 

u{0,y,z,t) = Unit) = (9/4- u(t)) • 16 • y{dy-y)/dl ■ z{d^-z)/dl ■ (1,0,0), 

with u{t) the average inflow velocity which is time-dependent, u{t) = sin(7rt/8). The considered time interval 
is [0, 8s], the viscosity is set to v = 10“^ s.t. the Reynolds number varies between Re = 0 and Re = 100. 


4-5.2. Drag and Lift. 

Again, the quantities of interest in this example are the (maximal and minimal) drag and lift forces Cd, 
cl that act on the disc. These are defined as 


with u 


max 




• Cx ds 



maX(g[o_8] = 1, r = 0.05, h = 0.41 and Fo the surface of the obstacle. 


Cy ds 
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Figure 5: Used mesh (left) and solution at t 


0.4 (right) for the benchmark problem in section [4.5| 


4-5.3. Numerical setup. 

We use an unstructured tetrahedral mesh consisting of 5922 elements, cf. Figure For the time 
discretization we use the modified fractional-step-0-scheme discussed in section |3.3| To compensate for the 
increase in the spatial accuracy for increasing polynomial degree fc, we adapt the number of time steps 
accordingly. The computations were carried out on a shared-memory computer with 24 cores. We comment 
on details of the computing times in remark El 



#ndof [K] 

max CD 

max Cl 

min Cl 

At [s] 

comp, time 

[s] 

k = 2 

169 

3.43046 

0.00262 

-0.016289 

0.0080 

492 X 

24 

fc = 3 

343 

3.29331 

0.00277 

-0.011099 

0.0080 

964 X 

24 

k — 4 

595 

3.29853 

0.00278 

-0.010762 

0.0040 

3 087 X 

24 

fc = 5 

939 

3.29798 

0.00278 

-0.011054 

0.0040 

6 670 X 

24 

ref. [59] 

11 432 

3.2963 

0.0028 

-0.010992 

0.01 

35 550 X 

24 

89 760 

3.2978 

0.0028 

-0.010999 

0.005 

214 473 X 

48 

ref. [SD] 

7 036 

3.2968 


-0.011 




ref. [58] 


[3.2,3.3] 

[0.002,0.004] 






Table 5: Numerical results for the benchmark problem “3D3Z” in m- Results obtained with different polynomial degrees and 
reference values. 


4 . 5 . 4 . Numerical results. 

In Table the results obtained are compared with the literature in terms of accuracy and computing 
time. We observe that we can achieve the same level of accuracy as the results in (and [SD]) with a 
computing time which is dramatically smaller. We note that the computing time per degree of freedom is 
actually worse than in |59j . Nevertheless, the same accuracy is achieved with much less degrees of freedoms 
using the high order method (fc > 2), s.t. our computations exceed the performance results in the literature. 
In the study [59] one of the conclusions is that their third order method {Q 2 /is much more efficient 
compared to lower order methods. We extend this conclusion in the sense that the use of even higher order 
discretizations, i.e. k > 2, increases efficiency even further. Moreover, high order discretizations can be 
implemented efficiently. One important component for the efficient handling of the Navier-Stokes equations 
with our high order discretization is the time integration using operator-splitting. 


Remark 11 (Computation times). We remark on the computation for the case k = 5. In that computa¬ 
tion approximately 65% of the computing time has been spend on the solution of linear systems for Stokes-type 
problems, 30% on convection operator evaluations and 5% on the setup and remaining operations. To solve 
the implicit convection problem (Step 2 in {4:1)) an average of 20 iterations has been applied. 


Remark 12 (Linear systems). In the test cases in this section we only applied direct solvers which is 
possible due to the (comparably) small size of the arising linear systems. For problems with increasing 
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complexity efficient linear solvers are mandatory. The development of suitable preconditioners of the Stokes 
problem is based on efficient preconditioning of the bilinear form Ah- For the scalar problem we could 
show poly-logarithmic bounds in k for the condition number of standard p-version domain decomposition 
preconditioners, cf. We plan to investigate suitable generalizations of this preconditioner for the Stokes 

problem in the future. 

5. Conclusion 

We presented and discussed a combined DG/HDG discretization tailored for efficiency. We summarize 
the core components. We split the Navier-Stokes problem into linear Stokes-type problems and hyperbolic 
transport problems by means of operator-splitting time integration. For the Stokes-type problem we use an 
i/(div)-conforming Hybrid DG formulation with a new modification: the projected jumps formulation. The 
Hybrid DG formulation facilitates the efficient solution of linear systems compared to other DG methods. 
The projected jumps formulation improves its efficiency even more. For the hyperbolic transport problem we 
apply a standard DG formulation. In numerical test cases we demonstrated the performance of the method. 
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